"""
电池生产质量控制分析
场景：某电池厂测试了两种生产工艺（工艺A和工艺B）生产的电池寿命
质量工程师需要验证：
1. 两种工艺的平均寿命是否有显著差异（均值检验）
2. 两种工艺的寿命稳定性是否相同（方差检验）
"""

import numpy as np
from scipy import stats
import matplotlib.pyplot as plt
import seaborn as sns

# 设置随机种子保证结果可复现
np.random.seed(2023)

# 模拟电池寿命测试数据（单位：小时）
# 工艺A：平均寿命=1200小时，标准差=50小时，样本量=50
process_a = np.random.normal(loc=1200, scale=50, size=50)

# 工艺B：平均寿命=1220小时，标准差=80小时，样本量=50
process_b = np.random.normal(loc=1220, scale=80, size=50)

print("===== 数据摘要 =====")
print(f"工艺A: 平均寿命={np.mean(process_a):.2f}小时, 标准差={np.std(process_a):.2f}小时")
print(f"工艺B: 平均寿命={np.mean(process_b):.2f}小时, 标准差={np.std(process_b):.2f}小时")

# ----------------------------
# 第一步：数据可视化
# ----------------------------
plt.figure(figsize=(12, 6))

# 箱线图比较  直观显示工艺B的中位数更高但离散程度更大
plt.subplot(1, 2, 1)
sns.boxplot(data=[process_a, process_b])
plt.xticks([0, 1], ['工艺A', '工艺B'])
plt.title('电池寿命分布比较')
plt.ylabel('寿命(小时)')

# 直方图比较  展示工艺B的寿命分布更分散，存在更多低寿命电池
plt.subplot(1, 2, 2)
sns.histplot(process_a, color='blue', alpha=0.5, label='工艺A', kde=True)
sns.histplot(process_b, color='red', alpha=0.5, label='工艺B', kde=True)
plt.title('电池寿命分布直方图')
plt.xlabel('寿命(小时)')
plt.legend()

plt.tight_layout()
plt.savefig('battery_life_comparison.png', dpi=300)
plt.show()

# ----------------------------
# 第二步：正态性检验（t检验的前提）
# ----------------------------
print("\n===== 正态性检验 =====")
# Shapiro-Wilk正态性检验
_, p_a = stats.shapiro(process_a)
_, p_b = stats.shapiro(process_b)

alpha = 0.05  # 显著性水平

print(f"工艺A Shapiro-Wilk p值: {p_a:.4f}")
print(f"工艺B Shapiro-Wilk p值: {p_b:.4f}")

if p_a > alpha and p_b > alpha:
    print("结论: 两组数据均符合正态分布，可以继续t检验")
else:
    print("警告: 数据不符合正态分布，建议使用非参数检验")
    # 此处为简化流程，我们继续使用参数检验，但实际生产应转为非参数方法

# ----------------------------
# 第三步：方差齐性检验
# ----------------------------
print("\n===== 方差齐性检验 =====")
print("目的: 检验两种工艺的电池寿命波动程度是否相同")

# 方法1: Levene检验（对非正态数据更稳健）
levene_stat, levene_p = stats.levene(process_a, process_b)
print(f"Levene检验统计量: {levene_stat:.4f}, p值: {levene_p:.4f}")

# 方法2: F检验（传统方法，但对正态性要求更严格）
f_value = np.var(process_b, ddof=1) / np.var(process_a, ddof=1)  # 计算F统计量
df1 = len(process_b) - 1  # 分子自由度
df2 = len(process_a) - 1  # 分母自由度
f_p = 2 * min(stats.f.cdf(f_value, df1, df2), 1 - stats.f.cdf(f_value, df1, df2))  # 双尾p值
print(f"F检验统计量: {f_value:.4f}, p值: {f_p:.4f}")

# 使用Levene检验结果做决策（更稳健）
if levene_p > alpha:
    print(f"结论(p={levene_p:.4f} > {alpha}): 两种工艺寿命稳定性无显著差异")
    equal_var = True  # 标记方差相等
else:
    print(f"结论(p={levene_p:.4f} <= {alpha}): 两种工艺寿命稳定性存在显著差异")
    equal_var = False  # 标记方差不相等

# ----------------------------
# 第四步：均值差异检验
# ----------------------------
print("\n===== 均值差异检验 =====")
print("目的: 检验两种工艺的平均电池寿命是否有显著差异")

# 执行独立样本t检验
# 根据方差齐性检验结果设置equal_var参数
t_stat, t_p = stats.ttest_ind(process_a, process_b, equal_var=equal_var)

# 计算置信区间
diff_mean = np.mean(process_b) - np.mean(process_a)
n1, n2 = len(process_a), len(process_b)
if equal_var:
    # 合并标准差
    pooled_std = np.sqrt(((n1-1)*np.var(process_a, ddof=1) + (n2-1)*np.var(process_b, ddof=1)) / (n1+n2-2))
    std_error = pooled_std * np.sqrt(1/n1 + 1/n2)
else:
    # 分别计算标准误
    std_error = np.sqrt(np.var(process_a, ddof=1)/n1 + np.var(process_b, ddof=1)/n2)

# 95%置信区间 (使用t分布的临界值)
dof = n1 + n2 - 2 if equal_var else (std_error**4) / ( (np.var(process_a, ddof=1)**2/(n1**2*(n1-1))) + (np.var(process_b, ddof=1)**2/(n2**2*(n2-1))) )
t_critical = stats.t.ppf(1 - alpha/2, dof)
ci_low = diff_mean - t_critical * std_error
ci_high = diff_mean + t_critical * std_error

print(f"t统计量: {t_stat:.4f}, p值: {t_p:.4f}")
print(f"工艺B与工艺A的寿命差异: {diff_mean:.2f}小时")
print(f"差异的95%置信区间: ({ci_low:.2f}, {ci_high:.2f})小时")

# 结果解读
if t_p > alpha:
    print(f"结论(p={t_p:.4f} > {alpha}): 两种工艺的平均寿命无显著差异")
else:
    print(f"结论(p={t_p:.4f} <= {alpha}): 两种工艺的平均寿命存在显著差异")

# ----------------------------
# 第五步：生产决策分析
# ----------------------------
print("\n===== 生产决策建议 =====")
print(f"工艺A寿命: {np.mean(process_a):.1f} ± {np.std(process_a):.1f}小时")
print(f"工艺B寿命: {np.mean(process_b):.1f} ± {np.std(process_b):.1f}小时")

# 分析质量稳定性
if equal_var:
    print("→ 两种工艺的寿命稳定性相同")
else:
    # 计算变异系数(CV)比较相对波动性
    cv_a = (np.std(process_a) / np.mean(process_a)) * 100
    cv_b = (np.std(process_b) / np.mean(process_b)) * 100
    print(f"→ 工艺B的波动性更大(CV: A={cv_a:.1f}% vs B={cv_b:.1f}%)")

# 分析平均寿命差异
if t_p > alpha:
    print("→ 两种工艺的平均寿命无明显差异")
else:
    if diff_mean > 0:
        print(f"→ 工艺B的平均寿命显著高于工艺A约{diff_mean:.1f}小时")
    else:
        print(f"→ 工艺A的平均寿命显著高于工艺B约{-diff_mean:.1f}小时")

# 综合生产建议
if not equal_var and t_p <= alpha:
    if diff_mean > 0:
        print("\n★ 推荐决策: 采用工艺B（寿命更长），但需改进工艺稳定性")
    else:
        print("\n★ 推荐决策: 坚持使用工艺A（更稳定），但可优化配方提高寿命")
elif equal_var and t_p <= alpha and diff_mean > 0:
    print("\n★ 推荐决策: 全面切换至工艺B（寿命更长且稳定性相同）")
elif equal_var and t_p > alpha:
    print("\n★ 推荐决策: 两种工艺效果相当，选择成本更低的工艺")
else:
    print("\n★ 推荐决策: 需要进一步实验优化工艺")

# ----------------------------
# 第六步：效应量计算（补充分析）
# ----------------------------
print("\n===== 效应量分析 =====")
# Cohen's d (标准化均值差异)
pooled_std = np.sqrt(((n1-1)*np.var(process_a, ddof=1) + (n2-1)*np.var(process_b, ddof=1)) / (n1+n2-2))
cohen_d = diff_mean / pooled_std
print(f"Cohen's d (效应量): {cohen_d:.2f}")

# 效应量解释
if abs(cohen_d) < 0.2:
    print("效应量: 微小 (差异在实际中可能不重要)")
elif abs(cohen_d) < 0.5:
    print("效应量: 小 (差异值得注意)")
elif abs(cohen_d) < 0.8:
    print("效应量: 中 (差异显著且有实际意义)")
else:
    print("效应量: 大 (差异非常明显)")

# 方差差异效应量
sd_ratio = np.std(process_b) / np.std(process_a)
print(f"标准差比率: {sd_ratio:.2f}")